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Abstract 



For an arbitrarily strong, spherically symmetric super-horizon curvature perturbation, we present 
analytic solutions of the Einstein equations in terms of the asymptotic expansion over the ratio 
of the Hubble radius to the length-scale of the curvature perturbation to set initial conditions for 
numerical computations of primordial black hole formation. To obtain this solution we develop a 
recursive method of quasi-linearization which reduces the problem to a system of coupled ordinary 
$L| ■ differential equations for the iVth order terms in the asymptotic expansion with sources consisting 

of a non-linear combination of the lower order terms. 
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I. INTRODUCTION 



The idea that large-amplitude matter overdensities in the Universe could have collapsed through 
self-gravity to form primordial black holes (PBHs) was first put forward by Zel'dovich and Novikov 
[H, and then independently by Hawking [4], more than three decades ago. This theory suggests 
that large-amplitude inhomogeneities in the very early universe overcome internal pressure forces 
and collapse to form black holes. A lower threshold for the amplitude of such inhomogeneities was 
first provided by Carr 0, 0] for a radiation-dominated epoch. The PBH contribution to the energy 
density increases with time during this epoch. For this reason, the PBHs formed considerably 
before the end of radiation-domination, possibly even before radiation-domination 0, aflect various 
cosmological and astrophysical processes even if their initial abundance is tiny. PBHs with mass 
smaller than ~ 10 15 g would have evaporated through Hawking radiation 0] and their abundances 
are constrained by big-bang nucleosynthesis 0413] and the gamma-ray background 13-15j], while 
holes with larger masses are constrained by dynamical and lensing effects [161 ] and by the stochastic 



gravitational wave background 17], [18|]. All these constraints are updated and summarized in 191 ]. 



Since the probability of PBH formation depends crucially on the statistical characteristics of 
the random field of primordial perturbations, PBHs provide a useful and unique tool to obtain 
independent constrains on the primordial power spectrum of inhomogeneities on extremely small 
scales which cannot be probed by any other methods. To make this cosmological tool more reliable, 
we must improve the prescription of the initial conditions 20|. Self-consistent initial conditions 
are very important for calculations of the probability of PBH formation [2l| and for relativistic 



hydrodynamical computations 2214251] . According to such computations the pressure gradients in 



the collapsing configuration play an extremely important role and are directly determined by the 
curvature profile in the initial configuration. 

Since PBHs can form only from highly non- linear curvature perturbations, the initial conditions 
of their formation must be consistent with the underlying non-linear theory such as the general 
relativity. The main objective of the present paper is to exclude the possibility that even highly 
sophisticated computer simulations could produce irrelevant results due to inconsistency of the 
initial conditions. For example, if we assume that the Universe outside the configuration is spatially 
flat, the mass of a perturbed configuration of radius r should be equal to the mass of unperturbed 
sphere of the same radius. As a result, a self-consistent density profile should be non-monotonic, 
i.e. along with the region of density excess, it should contain a region of density deficit, which 
drastically changes the effect of pressure gradients [22 ] . 

For an arbitrarily strong spherically symmetric super-horizon curvature perturbation, we present 
an analytic solution of the Einstein equations in terms of an asymptotic expansion over the ratio 
of the Hubble radius to the length-scale of the curvature perturbation under consideration. This 
method is similar to the gradient expansion 26-31] (previously known as the anti-Newtonian 
expansion [32]]) in spirit, but thanks to the spherical symmetry we can construct a solution to 
arbitrary higher order in this ratio from a single function characterisin g th e curvature profile. Note 



that the lowest-order solution in this program has been obtained in 20, l22l. 1331] . Here we develop the 
recursive method of quasi-linearization which reduces the problem to a system of coupled ordinary 
differential equations for terms of nth order in the asymptotic expansion of the density, pressure, 
velocity and metric, with sources which contain a non-linear combination of the lower order terms. 
Using this method, we obtain analytic expressions for all terms. 

The curvature profile, K\{r), to be defined below, appears in the source terms on the right- 
hand side of the relevant equations and deviations from homogeneity in density and velocity are 
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generated by the inhomogeneity of the curvature. To avoid confusion, we should say that these 
deviations do not involve small cosmological perturbations at all. Statistical characteristics of 
small perturbations are only relevant in this context when one calculates the probability of finding 
a configuration with a high amplitude perturbation of the metric. 

Dropping all terms of order greater than N, we obtain truncated asymptotic solutions of iVth 
order. Then for the arbitrary precision required by the intended accuracy and stability of the 
computer code, we obtain an upper limit on the time when such an iVth order truncated expansion 
can be used to set the initial conditions of fully non-linear numerical simulations of PBH formation. 
Later initial times obviously correspond to shorter computer runs. Thus our analytic solution helps 
to optimize numerical computations. 

The rest of the paper is organized as follows. In §11 basic equations are derived and in §111 
expansion coefficients are defined and their properties are described. Then in §IV we derive the 
recursive formulae for the coefficients in our problem and they are solved in §V for several specific 
initial curvature profiles. §VI is devoted to discussion and drawing conclusions. 



II. MATHEMATICAL FORMULATION OF THE PROBLEM 
A. The Misner-Sharp equations 

Assuming spherical symmetry, it is convenient to divide the collapsing matter into a system of 
concentric spherical shells and to label each shell with a Lagrangian comoving radial coordinate r. 



Then the metric can be written in the form used by Misner and Sharp 341 ] : 

ds 2 = -a 2 dt 2 + b 2 dr 2 + R 2 (d9 2 + sin 2 8d0 2 ), (1) 

where R, a and b are functions of r and the time coordinate t. We consider a perfect fluid with 
energy density p(r,t) and pressure p(r,t) and constant equation-of-state parameter 7, p(r,t) = 
7p(r, t). Expressing the proper time derivative of R as 

R /„\ 

« = -, (2) 

with a dot denoting a derivative with respect to t, we derive equations of motion for these variables 
as follows. 

First, from the (J!) component of the Einstein equations, we find 

b _ au' 

b~ R" [6) 

while the Euler equation yields 

°L = i—i. ( 4 ) 

1+7 p [) 

where a prime denotes differentiation with respect to r. We define the mass within the shell of 
proper radius R by 

fR(r,t) 

M(r, t)=4n p(r,t)R 2 dR, (5) 

J 

which gives 

M' = ^ P R 2 R'. (6) 



3 



Using ([3]), the ((j) component of the Einstein equations becomes 



R' 2 1 2 2GM 



and (J5J) can then be expressed as 



2 2GM\3 



M(r, t) = J p M + u* - — dV, (8) 

where dV = 4:7rR 2 bdr is the proper volume element. Equation (|8|) shows that M includes contri- 
butions from both the kinetic energy and the gravitational potential energy. Finally, combining 
©~C3), the evolution equation of M becomes 

M = -AirpR 2 R. (9) 
B. Quasi- homogenous asymptotic equations in new variables 

We consider the evolution of a perturbed region described by the above equations embedded in 
a flat Friedmann-Lemaitre-Robertson- Walker (FLRW) Universe with metric 

ds 2 = -dt 2 + S 2 (t)(dr 2 + r 2 d6 2 + r 2 sin 2 (10) 

which is a particular case of JT]). The scale factor in this background evolves as 

where t\ is some reference time. 

We denote the background solution with a suffix 0. In terms of the metric variables defined in 
([I]), we find 

a = 1, b = S(t), R = rS(t). (12) 



The background Hubble parameter is 



- ^ - 1 - f ■ < I3 > 



and the energy density is calculated from the Friedmann equation, 



We introduce a variable H defined by 



/ \ 3a 2 ,„ 
= 8rfF' (14) 



and another new variable H by rewriting H as 

H(t,r) = H (t)H(t,r). (16) 

The tilde-variable H, as well as most of the other tilde-variables introduced below, represent 
deviations of the solutions from the corresponding ones in the flat FLRW universe. Specifically we 
find 

a(t, r) = a (t)a(t, r) = d(t, r), (17) 
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b(t,r) = b (t)~b(t,r) = S(t)b(t,r), (18) 
R(t, r) = R (t)R(t, r) = rS{t)R(t, r), (19) 
p(t,r)=p (t)p(t,r)cxt- 4a p. (20) 

We define another variable Jx by 

M = —p R 3 p, (21) 
and the curvature profile K{t, r) is defined by rewriting b as 

K(t, r) vanishes outside the perturbed region so that the solution asymptotically approaches the 
background FLRW solution at spatial infinity. 

We denote the comoving radius of a perturbed region by n, whose precise definition will be 
given later, and define a dimensionless parameter e in terms of the square ratio of the Hubble 
radius Hq 1 to the physical length scale of the configuration, 

es Ufe) = =v " s2(i -" ) - (23) 

When we set the initial conditions for PBH formation, the size of the perturbed region is much 
larger than the Hubble horizon. This remains the case until the horizon mass becomes larger than 
the PBH mass. The horizon mass grows with cosmic time (for the radiation-dominated regime, 
the growth is directly proportional to time). This means e < 1 at the beginning, so it can serve 
as an expansion parameter to construct an analytic solution of the system ([2])- ([7]) to describe the 
dependence of all the above variables on the initial moment at which we set initial conditions. For 
the sake of brevity, below we will call this dependence "time evolution" . 

For our analytic expansion, it is convenient to rewrite the system of equations ([2])- ([7]) in terms 
of the tilde- variables, all of which tend to 1 at spatial infinity. From we find 

i + j^i = Map^)]' = 0, (24) 
a 1 + 7 p 

so 

a = F{t)p~^, (25) 

where F(t) is an arbitrary function of time. For convenience we choose F(t) = 1, then 

a = p"i+^. (26) 

Such a choice of F(t) corresponds to a frame of reference which is synchronous at spatial infinity. 

Since both p and a are positive definite, we may define p = In p and a = In a and then rewrite 
(USD as 

7 - 



1 + 7 

Using the tilde- variables, we can rewrite ([2]) as 



(27) 



so 



aR + tR = aaHR, (28) 



tR = aR($-l), (29) 
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where 

$ = aH. (30) 
Since R is positive, we can also define R by R = In R. Introducing a new time variable 

t = ta(£). (3D 

(f29|) is expressed as 

^=a(#-l). (32) 
From the definition of the curvature profile function K(r,t), (|22|) . we find 

Using ([2|) and ©, this can be rewritten as 

- r 2 ^ = 2#(1 - Kr 2 )^- = 2H(1 - Kr 2 )D r a, (34) 

(ri?)' 

where we have introduced the operator 

D r ^ (35) 
(rR)'dr 

We write the initial condition for (|34p as 

if(0,r) =Ki(r), (36) 

where K\(r) is an arbitrary function of r which vanishes outside the perturbed region, and define 
a new variable K by 

1 - K(t, r)r 2 = (1 - Ki(r)r 2 )K(t, r). (37) 

K is unity at spatial infinity like the other tilde- variables, but in contrast to the other tilde- 
variables, it describes the evolution of curvature deviation from the initial curvature profile rather 
than the deviation from the spatially flat Friedmann universe. Note that, from the definition (|22p 
of b(t,r), Ki(r) has to satisfy the condition 

Ki(r) < 1. (38) 

Physically, this condition ensures that the perturbed region does not form a closed universe which 
is causally disconnected from our universe 35, 36]. 

Differentiating ([37]) with respect to t and using (f34"|) . we find 



(1 - K- x r 2 )K = -r 2 K = 2H(1 - Kr 2 )D r a = 2H(l - K x r 2 )KD r ~a, (39) 

which yields 

tK = 2aHKD r a. (40) 



Since K is always positive, we can define another hat variable, K = InK, and rewrite (|40p using 
CUD as 

°K n f^n ~ 2aj$D r p 4 7 ~ , 

— = 2aHD r a = -—l JL = - ' $D r p. 41 

<9£ 1 + 7 p 3(1 + 
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Using (|2ip , we can write ([6]) in terms of the tilde- variables as 



47T 



— p Q R*~p ^ + 3 
3 \ n 



rR 



rR 



which leads to 



and hence 



fj, + 3/i — ~- = 3p — ~- 
rR rR 



P = P + g^r/i- 



In terms of the tilde- variables, ([9]) can be written as 



fl + p 



2 i? 

7 + 3 « 



p R 
37- „ = 0, 



so 



Defining 



6l) can be rewritten as 



tp, = 2p — 3aaH(p + 7/5). 
A + 7P 



/ 



1 + 7 ' 



In terms of H, the constraint equation ([7]) is expressed as 

8ttG „ ifr 2 



# 2 



"PoP 



and this gives 



:', '"' R 2 

eKr? 



p = H z + 



R? 



(42) 

(43) 
(44) 

(45) 
(46) 
(47) 

(48) 

(49) 
(50) 



One important property of the perturbation follows from (|49j) and the boundary conditions. 
The equation (|49p corresponds to the Priedmann equation of the flat FLRW universe 

8ttG 



Hi 



-Pa- 



using ([5]) and ([2"T]) . (j4*9|) and ([5T]) are combined to give 



H 2 ~ Hi 
Hi 



4tt 







P0 



R?Hf 



(51) 



(52) 



Noting the left-hand side and the second term of the right-hand side vanish at spatial infinity as a 
result of the boundary conditions and defining the energy density perturbation 



i(t,r) = ^ r) - w(tl = W , r) -i, 
Po(t) 

52]) leads to the following condition for 5: 

4vr / 5R 2 dR = 0. 



(53) 



(54) 

Namely, the mass excess in the center has to be compensated by the sorrounding mass deficit in 
order for the solution to coinside with the flat FLRW solution at spatial infinity. 

Equations (|2B]t, ([32"]) . ([IT]) . ([33]), ([35]) and ([50]) are the fundamental equations to solve. 
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III. EXPANSION OVER e AND QUASI-LINEARIZATION 

We now expand the tilde-variables over the parameter e as a first step to solving these funda- 
mental equations: 

oo 

X(t,r) = Y j e n (t)X in) (r). (55) 

n=0 

Note that by definition X(o) = 1 and XLs = 0. The hat-variables are expanded similarly. Differ- 
entiating this expansion with respect to t, one obtains 

oo „ oo 

X(t,r) = eY j ne n ~\t)X {n) {r) =^n £ «(t)I (fl) (r) ) (56) 

n=0 n=l 

hence 

;i) (n) = ¥*«,,. (57) 

Differentiating the expansion with respect to r, one finds 

oo 

*'(f,r) = X> B (r), (58) 

n=l 

hence 

{X') (n) =X{ n) . (59) 
Let us consider the product of two tilde- variables X\ and X 2 : 



OO 



Xl x 2 = £Vx 1(<) E' jjt m = EE^'W^ (60) 

\i=0 / \j=0 J i=0 j=0 

from which one finds 

n 

(XlX 2 )(n) = y^ y X l(i) X 2(n-i)- ( 61 ) 

Since -X"i(o) = -^2(0) = 1> one obtains 

(XiX 2 )( n ) = Xx( n ) + X 2 ( n ) + S( n )[X 1 X 2 ], (62) 

where 

rt-1 

S( n )[X 1 X 2 \ = y^Xx^^n-j). (63) 



Note that [X1X2] = S^[XiX 2 ] = 0. The most important feature of S^[XiX 2 ] is that it 
depends only on coefficients up to (n — l)th order. 

The relationship (|62p can be generalized to arbitrary functions of the tilde- variables. Let F\ 
and F 2 be arbitrary functions of tilde-variables and their time and space derivatives. Then 

n 

(F 1 F 2 ) {n) = *l(i)*2(n-i) = F l{n) F m + F m F 2(n) + S {n) [F^]. (64) 

i=0 

As a consequence of (f59|) and since X'^ = 0, one finds 

(X'F) (n) = X[ n) F {0) + S {n) [(X)'F], (65) 
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and 

(X 1 , X 2 / ) (n) = < S (r^) [(X 1 , )(X 2 , )]. (66) 
Similarly, using ([57]) and ([65]) and noting Z/ ) = 0) one obtains 

(XF) (n) = ^-X (n) F (0) + 5 (n) [(l)F] = ^ (x (n) F (0) + ^ ( ;,[IF]j , (67) 
where we have defined 

n-1 

S{ n) [XF] = mX m F {n _ m) . (68) 

m=l 

We also find 

(i l X2\ n) = - t Sl n) [X l X2]. (69) 

It is useful to obtain relationships between the expansion coefficients of the tilde-variables and 
those of the hat- variables. Suppose Y is some product of the positive-definite quantities a, R, p 
and K, such as 

Y = a pl R P2 p P3 K p \ (70) 
where pi, • • ■ ,P2 are integers. Then defining 

Y = lnY = p 1 a+p 2 R + P3P + PAK, (71) 
we can relate the expansion coefficients of 

1 d n Y 

Y(n) = -lKm— (72) 



and 



by 



y ' n\ e->o ae r ' 

1 , d n i 
'- (n) = — r hm — — 



1 d n Y 

Y (n) = -i^^r, (73) 



Y - V m! is f gm Y ^ Hm / <r " 

l / n ^— I'm. — 1 



' (m — 1)! e^o \m\de m J e->o \(n — m)lde r ' 
- ^ mY( m )Y(„_. m ). (74) 



n 

m=l 



Using Y( ) = 1, we obtain 

%)=%) + ^(n)[^]- (75) 



IV. EQUATIONS FOR ANALYTIC CALCULATIONS 

The fundamental equations to be solved in the following are ([57]). ([35]) . ([5T]) . ([35]) . ([35]) and 
()50p . We solve for the expansion coefficients X^ n ^ of each tilde or hat variable in the power series 
expansion with respect to e using these equations. To do this, we derive a set of recursive formulae 
to express Xir n \ in terms of Xjr m \ with m < n. 

First from (j27|) we find 

= ~Y^P(n), (76) 
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and 



which yields 

(1 + 7)a (n) + 7p (n) = ^S* {n) [p(p - a)]. (77) 

From (|48p and (|50p with (|77p . we can express -H"( n ) and /2( n ) in terms of the lower-order coefficients 
as follows. Using ([53jl and ([75]) . (j50 p leads to 

M(n) ~ 2 #( n ) = 5( n )[-H"-ff] + ^(e~ 2 ^)( n -i) 

+ rf (Ki - ^) {(e*) (B _D + (e- 2i V 1} + %-i)[eV 2 *]} 

= F (n) + W 1{n) (78) 

where 

F (n) = dirfKi - 2r 1 2 K iJ R (n _ 1) + rf ^ - 1^ , (79) 

W 1(n) = ^[Ml+r^fi-ijvDleV^] 

+ ^ (*i " ^) 5 (n-D ^ ~ 2r ? K l S (n-l) } • (80) 

On the other hand, (|48|) yields 

n/3A(n) = 2A(„) - 2($/) (n) . (81) 

Using the equalities 

*(n) = a(n) + H {n) + S( n )[5H], (82) 
A") = i^ilPW + £(»))> ( 83 ) 

we find 

($/)(n) = «(n) + ff(n) + 5( n )[a-H] + Yq--(M(n) +7P(n)) + £(n)[$/] 

= + rp/w + (i+ 7 )n 5 w^ " " )] 

+ S (n) [a£] + S (n) [*/|, (84) 
where we have used (|77|) in the last equality. From (|78|) and (|84|) we have 

hn) = ^ F{n) + Wl(n} ~ W ^ n )^ ^ 



H (n) 



1 



-— — — [A n (F {n) +W lin) ) + W 2{n) ], (86) 
where 

A n = Y^ f \( r r+l) n - r r\^ ( 87 ) 



W 2( n) = 2 (s (n) [aH] + S (n) m + n{1 \ 7) Sl n) \P(p ~ a)]) . (88) 

Now that we have expressed and #( n ) in terms of lower-order coefficients, we may use these 
coefficients to obtain recursive formulae for the other variables. For example, from ()44p we find 

T 

P(n) = A(n) + «A(n) + ^3(n), (89) 
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with 

W 3{n) = S (n) [(fi - ~ P ){tR)'] + T -S {n) \jl'R], 
where we have used (|65p . Then from (|77p we find 



T / T \ 

\Hn) + gM( n ) J + W4(n)> 



with 



Similarly ([32]) yields 



with 



From (HI 



with 



4(n) = "T^^H + (1 ^ ^J U W ~ *)]■ 



^ {n) ~ (l + 3 7 )n^ (n) + + W 5(n)), 



W" 5 („) = S {n) [aH). 



*M = - (l +7 )(l 7 + 3 7 )n (r/3 W + 
W 6(n) = rS (n) [/3'(*£)] + ( 2 + ^ + y) S*n)lK(rR)']. 



The corresponding tilde- variables are obtained from 



R(n) — R(n) + —S*n)[RR]> 

K{n)=K {n) + ±Sl ) [kK]. 

This completes our derivation of the recursive formulae. The first-order coefficients 
by the initial profile of the curvature inhomogeneity Ki(r) as 



r , w _3(l + 7)^iM 
/^(i)( r ) - ■ 



P(i)O) = 

R(i)(r) = - 



5 + 3 7 

r?Ki(r) 
~ 5 + 3 7 ' 

(l + 7 )rf(3ifj(r) + rif((r)) 
5 + 3 7 

7 rp(3Kj(r)+r^((r)) 
5 + 3 7 

(1 + 3 7 )rfl£i(r) + jrfrK[(r) 
5 + 18 7 + 9 7 2 ' 

2 7 r i 2 r(4Fq'(r) + riff (r)) 
(l + 3 7 )(5 + 3 7 ) 
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V. ANALYTIC SOLUTION 



Dropping all terms of order greater than N, we obtain truncated asymptotic solutions of Nth 
order. In order to use these solutions to set the initial conditions of numerical simulation of PBH 
formation, it is necessary to estimate an upper limit on time when such solutions are accurate 
enough to be used. Let the maximum acceptable error of the analytic solution be A. Then the 
latest epoch for which the analytic solution is accurate enough is determined by the first dropped 
terms of the truncated asymptotic expansions. If we use the asymptotic expansion of iVth order, 
the error of the asymptotic expansion for X is 

N 

£ 

n=0 



ERR £ e n X {n) )=X-J2 £ n X(n) = 0(e N+1 X (N+l) ). (105) 

\n=0 / 

Then we require 

e N+1 M (N+1) < A, (106) 

where M< n \ and M-y are defined by 



Mua = max 



and 



\m &i „Ms ,Mj> ,M Sf „Mb, „Mfi \ (107) 

= rnax|X (n) (r)|, (108) 

respectively. The error associated with the analytic calculation is less than A if it is calculated 
when 

e<e max = (109) 
y m (n+i) 

We solve the recursive relations obtained in §IV for four specific curvature profiles of the form 



Ki(r) 



1 + ^P 
2 Vo- 



exp 



1 fr 

2 W 



(110) 



where B describes slope of curvature profiles and a specifies the comoving length scale of curvature 
profile. Smaller values of B correspond to shallower profiles, and when B = the profile is simply 
Gaussian. The amplitude of the profile is set to unity at the origin where the same normalization 
is used as a spatially closed Friedmann universe in accordance with 20] . 

In order to represent the comoving length scale of the perturbed region, we use the comoving 
radius, n, of the overdense region. We can calculate r; by solving the following equation for the 
energy density perturbation defined by (f53j) : 

S(t, n ) = o. (in) 

Since the initial condition is taken at the superhorizon regime, when e is extremely small, the 
lowest-order solution ()10ip suffices to calculate n, which is obtained by solving 

m.in) + n K[{n) = o. (112) 

With the current choice of the functional form of K±{r), (|110|) . the solution of (|112|) is given by 

rf = 3a 2 for B = 0, (113) 
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and 



5jB _ 2 + ^(hB-2f + 2AB 
cj — for B^O. 



(114) 



We have obtained analytic solutions for curvature profiles with (B,a) = 
(1, 0.7), (0, 0.7), (1, 0.3), (0, 0.3), corresponding to wide and steep, wide and shallow, narrow 
and steep, narrow and shallow profiles, respectively. Plots of these profiles are shown in Figure [TJ 
Note that the physical length scale in the asymptotic Friedmann region is obtained by multiplying 
by the scale factor S(t), whose normalization we have not specified. We can therefore set up initial 
conditions for PBH formation with arbitrary mass scales by adjusting the normalization of S(t) 
which appears in the expansion parameter. 




FIG. 1: Initial curvature profiles K{(r) which are used as specific examples to obtain expansion 
coefficients of the tilde- variables. Note that these functions have to satisfy K\{r) < 1/r 2 . 

For these four specific profiles, expansion coefficients of the tilde- variables are calculated by 
solving the recurrence formulae ([85]) . ([86]) . (|89|) . (f9Tj) . (f9T|) and (j98|) numerically. Then, the quantity 
e max Was calculated for A = 10"\ 10~ 3 , 10 -5 and N = 1 - 7. The values of 

e max are summarized 

in Table HI When an asymptotic expansion of higher-order is used, e max is larger, so the analytic 
solution constructed is sufficiently accurate until a later time. For instance, one can see from the 
table that when an asymptotic expansion of first order is used, the numerical calculation has to 
be started at e = 0.0064 in order to maintain the accuracy of order 10 -5 , for the profile with 
(B,a) = (0,0.7). On the other hand, if we use an asymptotic expansion of seventh order, we 
can follow the evolution of perturbation until e = 0.51, maintaining the accuracy of 10 -5 , for the 
profile with (B,a) = (0,0.7). The dependence of e max on the order of the asymptotic expansion, 
N, is more clearly seen from Figure [21 in which the profile with (B,a) = (0,0.7) is used and A 
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(B,a) 


\ N 





1 


2 


3 


4 


5 


6 


7 


(1,0.7) 


1 n— 1 

10 


0.070 


0.30 


0.64 


0.73 


0.84 


0.93 


0.98 


1.0 


10 -3 


7.0 x 10~ 4 


0.030 


0.14 


0.23 


0.34 


0.43 


0.51 


0.57 


10 


7.0 x 10 


0.0030 


0.030 


0.073 


n i o 

0.13 


0.20 


0.26 


0.32 


(0,0.7) 


i n — 1 

10 1 


0.10 


0.64 


0.95 


1.3 


1.4 


1.5 


1.6 


1.6 


10 -3 


0.0010 


0.064 


0.21 


0.41 


0.56 


0.70 


0.82 


0.91 


10 


1.0 x 10 


0.0064 


0.044 


n i o 

0.13 


0.22 


0.32 


0.42 


0.51 


(1,0.3) 


1 n — 1 
10 


0.38 


1.4 


1.1 


1.1 


1.3 


1.5 


1.6 


1.7 


in — 3 

10 


0.0038 


0.14 


0.23 


0.36 


0.53 


0.70 


0.84 


0.95 


10~ 5 


3.8 x 10~ 5 


0.014 


0.049 


0.11 


0.21 


0.33 


0.44 


0.53 


(0,0.3) 


lO" 1 


0.56 


1.2 


1.6 


2.1 


2.5 


2.8 


3.2 


3.0 


10~ 3 


0.0056 


0.12 


0.34 


0.66 


1.0 


1.3 


1.6 


1.7 


10" 5 


5.6 x 10~ 5 


0.012 


0.074 


0.21 


0.40 


0.60 


0.85 


0.94 



TABLE I: Values of 

£max that satisfies the required accuracy of A for each pair of (B,<r) and 
different orders of asymptotic expansion, N. 

is set to be 10 _1 , 10~ 3 and 10~ 5 . The time dependence of A for asymptotic expansion of seventh 
order is shown in Figure [3l Note that it is determined by the first dropped eighth order terms of 
the expansions in this case. When the initial curvature fluctuation is wider and its profile steeper, 
expansion coefficients tend to be larger, so the errors in the analytic solution are also larger. 

Comparison of the time dependence of the errors associated with analytic solutions with different 
orders is shown in Figured The profile with (B, a) = (1, 0.7) was used for these plots. One can see 
clearly that the errors with higher-order expansions are relatively small and increase more slowly 
than those with lower-order expansions. Plots of the tilde-variables at e = 0.9, calculated using 
the asymptotic expansion of seventh order, is shown in Figure Note that errors associated with 
these plots are less than 10~ 3 from Figure [3l 

VI. DISCUSSION AND CONCLUSION 

In the present paper we have formulated a recursive method of quasi-linearization which can 
yield appropriate initial condition for PBH formation consistent with general relativity. The evo- 
lution of the profiles of the energy density perturbation 5 are shown in Figure [6l These profiles 
are calculated at e = 0.1,0.5 and 0.9. One can see that the region with 5 > 0, which corresponds 
to the central overdense region, is surrounded by the underdense region with S < so that (f54"|) is 
satisfied. 

We also introduce the averaged overdensity, denoted by 5 and defined as the energy density 
perturbation averaged over the overdense region as follows: 

(a \- x rR(t,r od ) 

S(t) = ( -nR(t,r od ) 3 \ J &tt5R dR. (115) 

Here r od (t) represents the comoving radius of the overdense region, which is numerically calculated 
from the solution of 5(t,r od ) = 0. It turns out that r Q d(t) is very close to n calculated from (j!12|) . 
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FIG. 3: The time dependence of errors associated with analytic solution obtained by asymptotic 
expansion of seventh order. 

i.e. lowest-order expansion. This feature can be directly observed in Figure 6, where the coordinate 
with 5 = hardly changes. 

The time evolution of the averaged overdensity 5 is shown in Figure [JJ For comparison, the 
results obtained using asymptotic expansions of first order are also shown. From the plots, one can 
confirm that higher-order corrections become more important as e gets closer to unity. When the 
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amplitude of initial curvature fluctuation is wider and its profile steeper, the density perturbation 
in the central region becomes larger, so that 5 tends to be larger. Therefore, it is more likely 
that wider and steeper initial curvature profiles lead to PBH formation after the perturbed region 
reenters the horizon. This confirms that considering the shape of profiles is crucial in the analysis 
of PBH formation. In addition, comparison of the time evolution of averaged overdensity 5 for 
N = 1 — 4 is shown in Figure [SJ The profile with (B,a) = (1,0.7) was used for these plots. 
When e <C 1, the plots coincide well with each other, but as e becomes larger, calculations using 
lower-order expansions start to deviate from those using higher-order ones. 

We have analyzed various configurations of curvature perturbations under the assumption of 
spherical symmetry to set up the initial condition for the numerical analysis of PBH formation in an 
optimal way with the help of the asymptotic expansion. In our analysis the curvature profile has a 
characteristic scale much larger than the Hubble radius initially, in accordance with the inflationary 
cosmology 37H39I] which predicts formation of superhorizon-scale curvature perturbations 40h43I|. 



This includes those perturbations which could lead to PBH formation 44145 71]. In a future paper 
we plan to calculate the probability of realization of the curvature profiles discussed above. Then 
we will eventually be able to relate the mass spectrum of PBHs with the parameters of inflationary 
models. 



Acknowledgments 

AGP acknowledges RESCEU for hospitality where this work was started. This work was sup- 
ported in part by JSPS Grant-in- Aid for Scientific Research No. 23340058 (JY), Grant-in- Aid for 



16 




FIG. 5: Profiles of the tilde- variables at e=0.9, which were calculated from an asymptotic expansion 
of seventh order, (a), (b), (c) and (d) were obtained from (B,a) = (1, 0.7), (0, 0.7), (1, 0.3), and 
(0, 0.3), respectively. Note that the errors associated with these profiles are less than of order 10~ 3 
from Figure [3l 
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FIG. 6: An illustration of "time evolution" of the density perturbation profiles. These were 
calculated when e = 0.1, 0.5 and 0.9. (a), (b), (c) and (d) were obtained from (B,a) = 
(1, 0.7), (0, 0.7), (1, 0.3), and (0,0.3), respectively. One can see that the region with 5 > 0, which 
corresponds to the central overdense region, is surrounded by the underdense region with 5 < 0. 
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FIG. 7: Time evolution of the averaged overdensity 5 for each of the initial curvature profiles. For 
comparison, the results obtained from the first-order asymptotic expansion are also shown by the 
dashed lines. 




N=4 
N=3 
N=2 
N=l 



0.0 0.2 0.4 0.6 0.8 



FIG. 8: Comparison of the time evolution of the averaged overdensity 5 for N = 1 — 4. The profile 
with (B,a) = (1,0.7) was used for these plots. 



19 



[1] Y. B. Zel'dovich and I. D. Novikov, Sov.Astron. 10, 602 (1967). 

[2] S. Hawking, Mon.Not.Roy.Astron.Soc. 152, 75 (1971). 

[3] B. J. Carr and S. Hawking, Mon.Not.Roy.Astron.Soc. 168, 399 (1974). 

[4] B. J. Carr, Astrophys.J. 201, 1 (1975). 

[5] M. Y. Khlopov and A. G. Polnarev, Physics Letters B 97, 383 (1980). 

[6] S. Hawking, Nature 248, 30 (1974). 

[7] Y. B. Zel'dovich, A. A. Starobinskii, M. Y. Khlopov, and V. M. Chechetkin, Sov. Astron. 
Lett. 3, 110 (1977). 

[8] I. D. Novikov, A. G. Polnarev, A. A. Starobinskii, and Y. B. Zel'dovich, Astron. Astrophys. 
80, 104 (1979). 

[9] B. V. Vainer and P. D. Naselskii, Astron. Zh. 55, 231 (1978), [Sov. Astron. 22, 138 (1978).]. 

[10] B. V. Vainer, O. V. Dryzhakova, and P. D. Naselskii, Pis ma Astronomicheskii Zhurnal 4, 344 

(1978), [Sov. Astron. Lett. 4, 185 (1978).]. 

[11] S. Miyama and K. Sato, Prog.Theor.Phys. 59, 1012 (1978). 

[12] K. Kohri and J. Yokoyama, Phys.Rev. D61, 023501 (2000), astro-ph/9908160. 

[13] D. N. Page and S. Hawking, Astrophys.J. 206, 1 (1976). 

[14] J. H. MacGibbon, Nature 329, 308 (1987). 

[15] J. H. MacGibbon and B. J. Carr, Astrophys. J. 371, 447 (1991). 

[16] B. Paczynski, Astrophys.J. 304, 1 (1986). 

[17] R. Saito and J. Yokoyama, Phys.Rev.Lett. 102, 161101 (2009), 0812.4339. 

[18] R. Saito and J. Yokoyama, Prog.Theor.Phys. 123, 867 (2010), 0912.5317. 

[19] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys.Rev. D81, 104019 (2010), 0912.5297. 

[20] A. G. Polnarev and I. Musco, Class. Quant. Grav. 24, 1405 (2007), gr-qc/0605122. 

[21] J. Hidalgo and A. Polnarev, Phys.Rev. D79, 044006 (2009), 0806.2752. 

[22] D. K. Nadezhin, I. D. Novikov, and A. G. Polnarev, Soviet Astronomy 22, 129 (1978). 

[23] I. D. Novikov and A. G. Polnarev, Soviet Astronomy 24, 147 (1980). 

[24] G. V. Bicknell and R. N. Henriksen, Astrophys. J. 232, 670 (1979). 

[25] M. Shibata and M. Sasaki, Phys.Rev. D60, 084002 (1999), gr-qc/9905064. 

[26] E. M. Lifshitz and I. M. Khalatnikov, Adv. Phys. 12, 185 (1963). 

[27] D. S. Salopek and J. M. Stewart, Class. Quantum. Grav. 9, 1943 (1992). 

[28] G. L. Comer, N. Deruelle, D. Langlois, and J. Parry, Phys. Rev. D 49, 2759 (1994). 

[29] Y. Nambu and A. Taruya, Class. Quantum. Grav. 13, 705 (1996). 

[30] I. M. Khalatnikov, A. Y. Kamenshchik, and A. A. Starobinsky, Class. Quantum. Grav. 19, 
3845 (2002). 

[31] Y. Tanaka and M. Sasaki, Prog. Theor. Phys. 117, 633 (2007). 

[32] K. Tomita, Prog. Theor. Phys. 54, 730 (1975). 

[33] I. Musco, J. C. Miller, and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009), 0811.1452. 

[34] C. W. Misner and D. H. Sharp, Phys.Rev. 136, B571 (1964). 

[35] B. J. Carr, T. Harada, and H. Maeda, arXiv:1003.3324 [gr-qc] (2010). 

[36] M. Kopp, S. Hofmann, and J. Weller, Phys. Rev. D 83, 124025 (2011). 

[37] K. Sato, Mon.Not.Roy.Astron.Soc. 195, 467 (1981). 

[38] A. H. Guth, Phys.Rev. D23, 347 (1981). 

[39] A. A. Starobinsky, Physics Letters B 91, 99 (1980). 



20 



[40] V. F. Mukhanov and G. Chibisov, Sov.Phys.JETP 56, 258 (1982). 

[41] A. H. Guth and S. Pi, Phys.Rev.Lett. 49, 1110 (1982). 

[42] S. Hawking, Phys.Lett. B115, 295 (1982), revised version. 

[43] A. A. Starobinsky, Phys.Lett. B117, 175 (1982). 

[44] J. Garcia-Bellido, A. Linde, and D. Wands, Phys. Rev. D 54, 6040 (1996). 

[45] H. M. Hodges and G. R. Blumenthal, Phys. Rev. D 42, 3329 (1990). 

[46] P. Ivanov, P. Naselsky, and I. Novikov, Phys. Rev. D 50, 7173 (1994). 

[47] J. Yokoyama, Astron. Astrophys. 673 (1997). 

[48] J. Yokoyama, Phys. Rev. D 58, 083510 (1998). 

[49] J. Yokoyama, Physics Reports 307, 133 (1998). 

[50] M. Kawasaki and T. Yanagida, Phys. Rev. D 59, 043512 (1999). 

[51] J. Yokoyama, Progress of Theoretical Physics Supplement 136, 338 (1999). 

[52] R. Saito, J. Yokoyama, and R. Nagata, Journal of Cosmology and Astroparticle Physics 2008, 
024 (2008). 

[53] A. Taruya, Phys. Rev. D 59, 103505 (1999). 

[54] B. A. Bassett and S. Tsujikawa, Phys. Rev. D 63, 123503 (2001). 

[55] A. M. Green and K. A. Malik, Phys. Rev. D 64, 021301 (2001). 

[56] M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama, Mod. Phys. Lett. A22, 1911 
(2007). 

[57] T. Kawaguchi, M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama, 
Mon.Not.Roy.Astron.Soc. 388, 1426 (2008), 0711.3886. 



21 



